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Axonal injury represents a critical target area for the prevention and treatment of trau- 
matic brain and spinal cord injuries. Finite element (FE) models of the head and/or brain 
are often used to predict brain injury caused by external mechanical loadings, such as 
explosive waves and direct impact. The accuracy of these numerical models depends on 
correctly determining the material properties and on the precise depiction of the tissues' 
microstructure (microscopic level). Moreover, since the axonal microstructure for specific 
regions of the brain white matter is locally oriented, the stress, and strain fields are highly 
anisotropic and axon orientation dependent. Additionally, mechanical strain has been iden- 
tified as the proximal cause of axonal injury, which further demonstrates the importance 
of this multi-scale relationship. In this study, our previously developed FE and kinematic 
axonal models are coupled and applied to a pseudo 3-dimensional representative volume 
element of central nervous system white matter to investigate the multi-scale mechani- 
cal behavior. An inverse FE procedure was developed to identify material parameters of 
spinal cord white matter by combining the results of uniaxial testing with FE modeling. 
A satisfactory balance between simulation and experiment was achieved via optimization 
by minimizing the squared error between the simulated and experimental force-stretch 
curve. The combination of experimental testing and FE analysis provides a useful analysis 
tool for soft biological tissues in general, and specifically enables evaluations of the axonal 
response to tissue-level loading and subsequent predictions of axonal damage. 
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INTRODUCTION 

An increase in the incidence and improvements in the diagnosis 
of traumatic brain injury (TBI) have increased awareness of the 
public about TBI's serious health effects. Axonal injury, which is 
considered to be a major contributor to cognitive dysfunction fol- 
lowing TBI, represents a critical focal area for TBI and spinal cord 
injury (SCI) prevention and treatment. FE models (FEM) of the 
head and/or brain are often used to predict brain injury caused by 
mechanical loading. Many attempts have been made to understand 
the injury mechanism and to define the mechanical parameters 
that govern axonal injury. Mechanical strain has been identified as 
the proximal cause of axonal injury (Wright and Ramesh, 2012), 
while secondary ischemic and excitotoxic insults associated with 
the primary trauma potentially exacerbate the structural and func- 
tional damage. Efforts to simulate the deformation, stress, and 
strain fields of brain tissues using computational models have pri- 
marily treated brain tissues as homogeneous materials, albeit with 
potential differences in properties between gray and white matter. 
As such, stress and strain fields were obtained in an average sense, 
which can be interpreted as macroscopic behavior (Kleiven, 2006; 
Ho and Kleiven, 2007; Cloots et al, 2008; Mao et al., 2010). 

There have been many reports on the mechanical characteriza- 
tion of central nervous system (CNS) tissue. For example, Miller 
et al. (2000) performed an in vivo indentation experiment to 



obtain a force-displacement curve, which was reproduced by using 
a three dimensional, non-linear FEM. By appropriately increas- 
ing material parameters describing the instantaneous stiffness of 
the tissue, they were able to reproduce the in vivo experimental 
curve. Shreiber et al. (2009) performed experiments to demon- 
strate the influence of myelin and glia on the tensile properties of 
the embryonic day 18 (E18) chick embryo spinal cord. They found 
that myelin and cellular coupling of axons via the glial matrix dic- 
tates the tensile response of the tissue. This work is an example of 
the influence of microstructural kinematics of CNS white matter 
that was introduced by Bain et al. (2003). The effects of non-linear 
material properties and the convoluted structure of the vascula- 
ture was investigated by Ho and Kleiven (2007) by using FEM 
and comparing the results with the 3D FEMs. Compared to a 
model without vasculature, the peak average strain reduction in 
a model with non-linear and linear elastic vasculature was 2 and 
5%, respectively, indicating that the influence of the vasculature 
on the dynamic response of the brain is minimal. 

Finite element (FE) analysis has become an important numer- 
ical tool for studying not only the mechanical behavior of brain 
tissues, but also prediction of TBI. Brain tissue has typically been 
characterized as a hyperelastic material (Meaney, 2003; Shreiber 
et al., 2009) and viscoelastic material (Horgan and Gilchrist, 2003; 
Post et al, 2012). For example, Post et al. (2012) used a linearly 
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viscoelastic material model combined with a large deformation 
theory to model brain tissue within a model of the head. The 
strain and Von Mises stress within the brain were extracted after 
simulating impacts that emulated typical loading curve shapes 
and peaks. However, the relationships between these bulk metric 
of the stress and strain fields within the brain relate to specific TBI 
neuropathologies and functional outcomes, such as concussion. 
Some efforts have been devoted toward this end. Sayed et al. pre- 
sented a biomechanical model for TBI and soft tissue damage. A 
variational constitutive model for soft biological tissues was uti- 
lized to reproduce axonal damage and cavitation injury through 
inelastic deformation (Sayed et al., 2008). Mechanical damage of 
brain tissue was classified as volumetric (compression/tension) 
and shear-type which were associated with intracranial pressure 
and shear stresses, respectively. 

The aforementioned FEM work is based on homogenized 
assumptions of brain tissue properties. However, the orientation 
and kinematic behavior of the microstructural constituents of 
white matter - including axons and the surrounding glia and 
vasculature - introduce significant anisotropy and complex bulk 
behavior. Thus, the accuracy, and perhaps ultimately the utility 
of FEM simulations depends not only on correctly determining 
the material properties, but also on the precise depiction of the 
tissues' microstructure (microscopic level including axon, glial 
cell, and blood vessels, etc.), to effectively relate the bulk tissue 
response to the actual elements that are injured. Recent compu- 
tational studies have recognized the potential importance of this 
multi-scale behavior. Karami et al. (2009) performed FE analysis 
(FEA) using a unit-cell comprised of undulated axons embedded 
within and fully coupled to a glial matrix. Their micromechanical 
FEM provides insight into how the local stress and strain fields 
are affected by the tortuosity (which we use interchangeably with 
undulation) of an axon, which is defined as the total axon path 
length divided by its end-to-end length. However, experimen- 
tal characterizations by Bain et al. (2003) and Hao and Shreiber 
(2007) clearly demonstrate that axons are not fully coupled to the 
glial matrix, and instead demonstrate complex, strain-dependent 
coupling behavior. Pan et al. (2011) developed a unit-cell model 
that included multiple axon segments with different tortuosities 
and coupling that emulated experimental data (Bain et al, 2003). 
The extended unit-cell model captured the kinematic behavior 
of axons very well when compared with experimental results. 
They found that the degree of coupling within the tissue affects 
the continuum mechanical properties of the tissue. Cloots et al. 
(2011) developed another micromechanical model to study how 
axonal strain is affected by the inclusion of other local cells. Var- 
ious factors, such as the inclusion of stiffness with respect to the 
surrounding tissue, the axonal orientation, and the maximum 
diversion angle of the axons, significantly influenced the axonal 
strain. 

Although the extended unit-cell model by Pan et al. (2011) 
successfully captured the kinematics behavior of axons, the model 
was simplified and did not effectively allow representation of the 
wide distribution of axon tortuosity that is evident in white mat- 
ter. To overcome this limitation, a pseudo-3D model is proposed 
in this study. A pseudo-3D representative volume element (RVE) 
is built using a random process. Its axons, resembling the axonal 



structure of CNS white matter, are embedded into an extracellu- 
lar matrix (ECM). The kinematic axonal coupling model, which 
demonstrated good correlation to experimental observations in 
our previous study, is applied to the RVEs in the FEM, and the 
tissue-level mechanical response is investigated. The FE model 
employs our previously published approach to introduce strain- 
dependent kinematic coupling of axons to matrix and published 
experimental results for mechanical properties. 

MATERIALS AND METHODS 

An inverse FE procedure is used to identify the material parame- 
ters of spinal cord white matter by combining uniaxial testing and 
FEM. An objective function is selected for minimization through 
which material parameters are determined. Hereafter, the methods 
employed to this end are presented. 

REPRESENTATIVE VOLUME ELEMENT OF CNS WHITE MATTER 
RVE construction 

The FE model of an RVE is developed using ABAQUS 6. 1 1 and 
Python scripting. An RVE is a sub-domain of a composite material 
whose size needs to be sufficiently large to include necessary fiber 
information such that the macroscopic homogeneous mechani- 
cal behavior of the composite material can be derived. To reduce 
computational demands in a FEA, one might need to limit the 
size of an RVE while retain necessary fibers. The interested reader 
is referred to Iorga et al. (2008). In this study, a pseudo-3D RVE 
for the CNS white matter is composed of many undulated axons 
that are embedded in the ECM. A cuboid (x = 0.4 |xm, y = 10 u,m, 
z = 5.68 |xm) is created representing the matrix surrounding the 
axons as shown in Figure 1. The axons are undulated and their 
geometry is complex. The undulation varies from axon-to-axon 
as well as along an individual axon. The distribution of undula- 
tion of the axons within the RVE model is based on those found by 
Bain et al. (2003). An undulated axon is represented by a poly-line 
and a cubic spline: a poly-line is composed of multiple connect- 
ing line segments extending from random points on one face of 
the cuboid (z = 0 plane) to its opposite face (z = 5.68 plane); and 
the cubic spline, representing the backbone of an axon, is gen- 
erated by depicting every key point on the poly-line. Sweeping a 
circle along the cubic spline then generates the axons. In this study, 
the diameter of an axon is fixed to 0.4 u,m, although it may vary 
from axon-to-axon in the range of sub-micron to about 20 |im for 
porcine optic and sciatic nerves (Assaf et al., 2008). The portion 
of an axon that falls outside of the cuboid is trimmed. The volume 
fraction of axons is set to 53% (Karami et al, 2009). The total 
number of axons in the matrix is 33, which is determined by the 
volume fraction, as shown in Figure 2. The average undulation of 
the axons ranges from 1.05 to 1.25. 

Kinematic of CNS white matter 

In this pseudo-3D RVE model, the axons interact with the sur- 
rounding matrix via the embedded element approach in ABAQUS. 
Embedded elements allow the introduction of an axon fiber within 
a matrix without removing the matrix material where the fiber 
resides. If the elements representing the fiber are not specified as 
embedded, there is no interaction between the matrix and fiber. 
In our model, this coupling is governed by the axon transitional 
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kinematic model (TKM) (Pan et al., 201 1), a model of axon-ECM 
interaction derived from statistical analysis of experimental data 
(Bain et al., 2003). Specifically, a single axon was divided into 50 
sub-segments along the axial direction. According to the TKM, the 
percentage of an axon that is embedded in the matrix depends on 
its current undulation and applied stretch (Pan et al., 2011). In 
general, there is a transition undulation that dictates whether an 
axon is coupled or not. As stretch is applied, the axons straighten 
and undulation decreases. Accordingly, the likelihood of reaching 
the transition undulation, thereby coupling the deformation of 
the axon to the matrix, increases with increasing stretch. If the 
undulation of an axon is less than 1.08, the percentages of that 
axon being coupled are 8, 20, and 44% of its whole length when 
the corresponding applied stretches are less than 1.06, 1.12, and 
1.25, respectively. If undulation is greater than 1.08, the axon does 
not interact with the matrix (except both ends are tied to the 
matrix at the RVE surfaces). The threshold undulation, 1.08, was 
determined from experimental data by Bain et al. (2003). 

To demonstrate our model, three RVEs (RVE1, RVE2, RVE3) 
with randomly generated axons were simulated. Initially, the mod- 
els were loaded with applied stretches ranging from 1.0 to 1.06 in 
the z-direction, which was done by fixing U z = 0 on the Z = 0 
face and applying the displacement, [7 za , on the opposite face. 
Next, this was done with a range of 1.06-1.12 and finally with a 
range of 1.12-1.25. The tortuosity of each axon at the initial stage 
(X = 1.0) and after each step (X = 1.06, 1.12, and 1.25) was com- 
puted to compare to the experimental data and previous modeling 
approaches. Interactions with embedded elements were updated 
after each step. 

INVERSE FE PROCEDURE 

An inverse FE procedure was developed to identify material para- 
meters of spinal cord white matter by combining uniaxial testing 
and FEA using the pseudo-3D RVE. The material properties of 
axons and ECMs were modeled with Ogden's isotropic, large defor- 
mation, hyperelastic material model in the strain energy form 
(Meaney, 2003): 

W-=^(X? + \« + X*-3) (1) 
or 

where \l represents the shear modulus, a is a material-dependent 
parameter that introduces non-linear behavior, and X; are the 
principal stretches. 



A satisfactory balance between simulation and tensile-testing 
based experimental data was achieved via an optimization process. 
The squared error between the simulated and experimental 
data was used as an objective function to be minimized. The 
stress-stretch curve obtained from an experiment is denoted as: 

^exp = ^exp (X) (2) 

where S exp is the engineering stress and X is the stretch of a sample 
from an uniaxial tensile test. Assuming the same experiment was 
simulated with FEM, the simulated stress-stretch curve would be 
denoted as: 

Ssim = S sim (>.,ai,...,o( n ) (3) 

where S s i m is the simulated stress, X is the prescribed stretch, and a\ 
through a n are the material parameters to be determined. To opti- 
mize the material parameters, an error function, y — S eX p — 5 s i m , 
is defined. The objective function to be minimized is y 2 at pre- 
established values of the independent variable, X, that fall within 
the test range. A practical approach is to discretize the stress-stretch 
curve and use a least-squared method to determine the parame- 
ters (Zhang and Gan, 201 1). To demonstrate the application of the 
Pseudo-3D RVE, a single-variable example regarding a quasi-static 
uniaxial tensile-testing based experimental stress-stretch curve 
(Shreiber et al., 2009) combined with the FE simulation was cre- 
ated to determine the shear modulus, \i, of axons. It was assumed 
that both the axons and ECM follow the Ogden hyperelastic model 
in Eq. (1) (Meaney, 2003). The non-linear parameter, a , was set 
to be the nominal value, 8.22, obtained in Shreiber et al. (2009) 
via regression analysis of the experimental data. For simplicity, the 
same value is used for both the matrix and the axons, although 
it might be different for the axons and the matrix in general. The 
summation of the squared error of the stress-stretch curve at dif- 
ferent pre-established values is used as the objective function. The 
problem can be described as: 

Min E (|X) = ^. = 1 (Sexp_, - S sim _,-) ^ 
S.t. flmin — \^ — H-max 

where the squared error, E([i), is the objective function, shear 
modulus, [i, is the material parameter to be optimized, S ex p_i and 
Ssimj are experimental and simulated stress at stretch level X 1 
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S, Mises 




FIGURE 2 | Numerical uniaxial tensile experiment at the stretch level of 1.06. (A) von Mises stress contour of the matrix of the RVE in which the undulated 
axons are embedded. (B) von Mises stress contour of the axons. Highly localized stress implies the effect of the axon geometry on local stress and strain field. 
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the i pre-established stretch. Note that this is different from the 
principal stretch Xj. 

The optimization for material properties is described in the 
flow chart in Figure 3. The optimization loop was achieved with 
a program coded in Python. The FE solver was used to solve the 
uniaxial tensile test numerically. The FE results were processed 
by the optimization program and compared to the experimental 
data. A new guess of parameters was generated based on the cur- 
rent comparison. The iteration is terminated when the optimized 
parameters were found using the tolerance of 1.0 kPa in the golden 
search algorithm. 

RESULTS 

The von Mises stress contours of the ECM and the axons subjected 
to a stretch rate of 1.06 in the z-axis are shown in Figures 2A,B, 
respectively. Since the interaction between matrixes and axons 
depends on the geometry of the axons, the stresses are highly local- 
ized. The overall average uniaxial stress-strain curves are plotted in 
Figure 4, where S33 is the nominal engineering stress component 
in the z-direction. 

The axons' tortuosity is plotted against the applied stretch in 
Figure 5. It is shown that the undulated axons become less wavy 
when the applied stretch increases. In this pseudo-3D RVE model, 
the tortuosity decreased faster at smaller applied stretches as com- 
pared to other simulation models. The experimental data from 
Bain et al. (2003) and the simulation data from Karami et al. 
(2009), who employed a unit-cell method with fully constrained 
coupling between axons and glia, are included for comparison, 
as well as our results from our previous TKM using a unit-cell 
approach (Pan et al., 2011). Clearly, fractional coupling between 
the axons and glia yields more accurate stretch-undulation data 
than Karami et al.'s model with perfectly bonded axons and 
glia. The current pseudo-3D RVE approach also gives consis- 
tent axon kinematics. Curves obtained from RVE1 and RVE3 
in this study are close to that from the unit-cell model and 
experimental measurements. The curve from RVE2 is slightly off 
the trend since its initial undulation was slightly higher than 
1.13. In this case, the evolution of axonal coupling to matrix 
with stretch is delayed. A drawback of these simulations is that 
the percentage of axon embedment was held constant and only 
updated at discrete stretch levels, whereas in reality axon kinemat- 
ics transition continuously. We expect that allowing the constraints 
to evolve with stretch will allow us to capture the non-affine, 
independent behavior of axons and glia at lower stretch levels 
and enable a more rapid change in undulation at low stretch 
levels. 

In order to determine the material properties of white matter in 
tension, an inverse FE simulation using the RVE and experimental 
data was performed. The controlled experimental stress-stretch 
data of chick spinal cords in Shreiber et al. (2009) were extracted 
and plotted in Figure 4 for inverse FE simulation. The non-linear 
parameter a in the Ogden hyperelastic model was curve-fit to be 
8.22 ± 1.27, which is used for both the axon and the ECM in the 
inverse simulation. The shear modulus of axon is assumed to be 
three times as stiff as that of ECM (Meaney, 2003). 

The stress-stretch curve was sampled at 20 equidistant intervals. 
The squared error was minimized with a golden search algorithm 
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FIGURE 3 I Flow chart of the optimization process 



for this univariate problem in order to find the shear modulus of 
an axon, [i. An initial guess for [i was given for the algorithm. A 
pre-determined solution bracket (20, 50 kPa) was also provided 
to shorten the computational time. The solution tolerance was 
set to 1.0 kPa for the unknown parameter. The procedure was 
automated using Python scripts. After nine iterations, the opti- 
mized parameter \l = 33.28 kPa was obtained using the nominal 
value a = 8.22. The residual of each iteration was tabulated in 
Table 1 and plotted in Figure 6. It is shown that the unknown 
parameter converges to the optimized value as the residual reaches 
its minimum. The stress-stretch curves obtained in all iterations 
were plotted in Figure 7. The simulated stress-stretch curves con- 
verge to the experimental curve. To investigate the sensitivity of 
the results to the parameter a, the optimization process was also 
performed at a = 6.95 and 6.49. The respective optimized para- 
meters are (i = 36.63 and 29.64 kPa and the respective residuals are 
tabulated in Table 1. The optimized parameters reside within the 
curve-fitting range of 32.8 ± 9.53 kPa obtained from the experi- 
mental data (Shreiber et al, 2009). It can be concluded that shear 
modulus is not sensitive to the non-linear parameter a since the 
change of \i is less than 11% and it is within the simulated range. 
The stress-stretch corresponding to the optimized parameter of 
|x = 33.28 kPa is also plotted in Figure 4 showing good agreement 
between the curves from the experiment and the optimization 
process. 
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FIGURE 4 | Stress-stretch curves of E18 chick embryo spinal cord: the converged simulation vs. experimental in Shreiber et al. (2009). 




FIGURE 5 | Predicted evolution of undulation with stretch. Simulation results are compared to experimental results from Bain et al. (2003) and to simulation 

results from Karami et al. (2009) and Pan et al. (2011). 
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DISCUSSION 

A novel pseudo-3D RVE with axon kinematics was developed to 
capture the kinematics of axons and the stress-stretch behavior of 
CNS white matter. The advantage of the embedded element model 
is that the complex geometry of the surrounding matrix is simpli- 
fied since the volume of axons was not excluded from the matrix. 
The kinematics of axons can be retained using the simplified axon 



Table 1 | Sensitivity analysis on the non-linear parameter a. 
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Change of the optimized parameter is within the range of ±11% when the para- 
meter a varies within the experimental error range. E(\l) is the squared error 
between the experimental and the simulation values. 
Bolded values are final results for \l with varying a. 



kinematic model, where a percentage of axons are transitioned 
from uncoupled to coupled based on the undulation and applied 
stretch. Stresses and strains are highly localized depending on the 
local undulation of axons and, importantly, on the interaction of 
the axon with the matrix. The effects of complex axon geometry 
on the local stress and strain fields are significant. The local stress 
and strain fields are highly heterogeneous nearby the axon, even 
though a uniform global strain field (£33 = 0.06, 0.12, and 0.25 in 
this case) is applied to the boundary of RVE. This implies that 
if a mechanical strain on axons is used as an injury criterion, a 
macroscopically homogeneous model is not sensitive enough to 
capture the local variation and a local "sub-model" or "micro- 
model" has to be considered and implemented. With the help of 
high resolution MRI techniques (McNab et al, 2009), one maybe 
able to obtain axon orientation of a human brain and build local- 
ized RVEs. The RVE will be useful to provide detailed, local axon 
deformation at the tissue-level. More work needs to be done to 
incorporate the model into a tissue-level representation of brain 
tissue white matter along this line. 

The application of a pseudo-3D RVE was demonstrated in 
the inverse optimization of the shear modulus of a hyperelastic 
model for the chick embryo spinal cord. Limited data is avail- 
able for the stress-strain properties of white matter in tension, 
and the chick embryo spinal cord studies clearly demonstrated 
that (1) axon kinematics transition from uncoupled to coupled 
with increasing stretch (Hao and Shreiber, 2007); and (2) cou- 
pling of axons to the glial matrix significantly affects the tensile 
stress-strain response (Shreiber et al., 2009). The E18 chick embryo 



300 



250 - 



200 



1 150 



100 





• 






■ - - ■ 




















★ 

■ 










opt 


mized va 


lues 








■ 


■ 


a 


= 6.95 


0 


0 


a 


= 8.22 






a 


= 9.49 



45 



50 



Optimization parameter 



FIGURE 6 I Plot of function value at various shear modulus. It shows that the squared error of the experimental curve and the simulated curve is minimized 
at |i = 36.63, 33.28, and 29.64 kPa for a = 6.95, 8.22, and 9.49, respectively. 
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FIGURE 7 | Inverse simulation example: stress-stretch curves of chick embryo spinal cord with a = 8.22 The arrow points to the stress-stretch curve 
corresponding to the optimized parameter |a = 33.28 kPa. 



spinal cords include gray and white matter, as well as the meningeal 
layers, and each may potentially contribute to the mechanical 
properties. If the spinal cord when tested in uniaxial tension is 
approximated as three concentric mechanical elements in a par- 
allel composite, the dura, arachnoid, and pial maters are too thin 
in the chick embryo to provide any significant contribution to the 
mechanical response. At El 8, the white matter comprises approx- 
imately 56% and the gray matter 44% of the cross-sectional area 
of the chick embryo spinal cord (Hao and Shreiber, 2007), so both 
tissues may contribute to the bulk tensile properties. However, 
Miller and Chiznei reported that brain tissue, with a predom- 
inance of gray matter with white matter that was not aligned 
with the direction of stretch, reaches a peak stress which is more 
than an order of magnitude lower than white matter in tension 
(Miller and Chinzei, 2002). Moreover, they found that stress- 
strain response of the brain tissue in tension was hypoelastic, 
contrary to the hyperelastic response of the chick embryo spinal 
cord (Shreiber et al, 2009) as well as the properties of the rat 
spinal cord (Fiford and Bilston, 2005). As such, the tensile prop- 
erties of the spinal cords are expected to be dominated by the 
response of the white matter. Unsurprisingly, the simulated stress- 
stretch responses correlate well with the experimental results. The 
parameter obtained through the inverse FE optimization proce- 
dure agrees well with the fitted material parameter. It is noted that 
while only a univariate problem is shown in this paper, a non- 
linear least-squared method can be used for multiple-parameter 
problems. 



The purpose of this study was to develop and validate a multi- 
scale computational approach that captures the microkinematic 
response of axons and the tissue-level response of the material 
in tension. The developing chick embryo spinal cord was cho- 
sen as our model tissue based on our previous characterizations. 
This chick embryo spinal cord conferred several advantages. Chick 
embryo tissue is readily available through procurement of fresh, 
fertilized eggs and is cost effective when compared to small mam- 
mals. Maintaining eggs during development is relatively simple, 
as is harvesting tissue at a desired developmental stage. Tissue can 
be rapidly harvested and tested - the experiments described in 
Shreiber et al. (2009) and Hao and Shreiber (2007) - were com- 
pleted within an hour on tissue harvest, well within the window 
described in Gefen et al. (2003) and Gefen and Margulies (2004) 
for correlation of ex vivo to in vivo properties of CNS tissue. Similar 
to the optic nerve, which was tissue used for the original descrip- 
tion of axon kinematics (Bain et al., 2003), axons in the spinal cord 
are predominantly oriented axially, the spinal cord also provides 
a system in which the prevailing direction of axons is parallel to 
the direction of stretch, which allows a direct interpretation of the 
relationship between macroscopic tissue stretch and microscopic 
axonal deformation. Along these lines, the fundamental behaviors 
observed in the chick embryo model paralleled those in the adult 
guinea pig optic nerve: undulated axons that straightened with 
stretch and demonstrated a transition from non-affine behavior 
at low stretch levels to affine behavior at high stretch levels. In 
addition, the form of the bulk stress-strain response of the chick 
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tissue in tension matched that from rat spinal cord (Fiford and 
Bilston, 2005) and the adult guinea pig optic nerve (Bain et al, 
1997; Bain and Meaney, 2000). As such, we believe that the model 
described herein is an accurate representation of CNS white mat- 
ter in tension that can be applied to multi-scale analysis of TBI 
and SCI. 

From the composite material mechanics of point of view, 
the mechanical properties of a composite heterogeneous mate- 
rial depend strongly on the volume fraction of the fiber; the 
higher the fiber volume fraction, the stiffer the composite mater- 
ial. In this study, the CNS white matter was treated as a biological 
fibrous composite material comprised of axons (fibers) and ECM 
(matrix). Changes of the average axon diameter will not affect the 
current results as long as the volume fraction remains unchanged. 
This is because the X- and Y-dimensions of the RVE need to be 
scaled proportionally to the axon diameter in order to maintain the 
same volume fraction. In this paper, we only consider a uniform 
axon diameter. It will be straightforward to adapt the methodol- 
ogy into the RVE with axon diameters that follow a distribution. 
In that case, a relative large number of axons are needed and the 
size of the RVE will increase. 

To fully implement this multi-scale approach for analysis of 
human TBI or SCI, a number of parameters should be iden- 
tified and/or tested. First, the tortuosity of axons in different 
regions of the brain and for different ages must be character- 
ized. The chick embryo and guinea pig studies demonstrated that 
tortuosity changes with development and can differ between tis- 
sues/locations and/or species. Second, the RVE approach should be 
validated against other modes of mechanical loading. We focused 
on the multi-scale response of axons in white matter to ten- 
sion because tensile strain has been recognized as the proximal 
cause of axonal damage (Margulies and Thibault, 1992). How- 
ever, tension can arrive from any multi-axial, deviatoric loading 
condition, and it is important that the RVE approach accurately 
emulate the bulk response of tissue in shear and compression. 
Finally, the RVE should be used to simulate a model of injury 
to identify and compare to existing injury criteria; similarly to 
the guinea pig optic nerve characterization study (Bain et al, 
2003) which was combined with in vivo injury studies (Bain and 
Meaney, 2000) to evaluate axon- and tissue-level injury criteria 
for functional damage, which was measured using electrophys- 
iology, and morphological damage, which was identified with 
immunohistochemistry. 
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